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ABSTRACT 

Numerical simulations of star formation frequently rely on the implementation of 
sink particles, (a) to avoid expending computational resource on the detailed internal 
physics of individual collapsing protostars, (b) to derive mass functions, binary statis- 
tics and clustering kinematics (and hence to make comparisons with observation), and 
(c) to model radiative and mechanical feedback; sink particles are also used in other 
contexts, for example to represent accreting black holes in galactic nuclei. We present 
a new algorithm for creating and evolving sink particles in SPH simulations, which 
appears to represent a significant improvement over existing algorithms - particularly 
in situations where sinks are introduced after the gas has become optically thick to its 
own cooling radiation and started to heat up by adiabatic compression, (i) It avoids 
spurious creation of sinks, (ii) It regulates the accretion of matter onto a sink so as 
to mitigate non-physical perturbations in the vicinity of the sink, (iii) Sinks accrete 
matter, but the associated angular momentum is transferred back to the surrounding 
medium. With the new algorithm - and modulo the need to invoke sufficient resolution 
to capture the physics preceding sink formation - the properties of sinks formed in 
simulations are essentially independent of the user-defined parameters of sink creation, 
or the number of SPH particles used. 

Key words: stars: formation - galaxies: nuclei - methods: numerical - hydrodynam- 
ics - gravitation 



1 INTRODUCTION 

Star formation is a co llective process: many stars are born in 
multiple systems (e.g. iDuquennov"^ Mayor 1991), and most 
are born in clusters (e.g. lLada fc Ladall2003h . Consequently, 
a realistic simulation that captures the important interac- 
tions between neighbouring protostars, as they form, is only 
feasible if some of the detailed internal physics of individual 
protostars is sacrificed. Otherwise all the available compu- 
tational resource is commandeered by the first protostar to 
form, and events in the rest of the computational domain 
grind to a halt. It is for this re ason that Lagrangj an sink 
particles were first developed bv lBate et al.l (| 19951 ). Subse- 
quently they have also been used to model single accreting 
objects, for example black holes in galactic nuclei, where the 
focus is not on the detailed hydrodynamics in the immediate 
vicinity of the black hole, but on feedback from th e black 
hole on much larger scales (e.g. ISpringel et al]|2005h . 

In the most basic implementation, a sink is created 
when and where some lump in the computational domain 
looks like it will inevitably condense out gravitationally (into 
a single star, or a tight multiple system, or a black hole). 



Sinks can also be introduced into the initial conditions for a 
simulation, to represent pre-existing condensations (stars or 
black holes). A basic sink is a point mass, and therefore it 
interacts with the surrounding matter gravitationally. Once 
created, a sink can grow by accreting nearby matter that 
looks like it would inevitably become part of the same con- 
densation. Even at this most basic level, there are issues 
with how to identify lumps that should be converted into 
sinks, and how to determine which matter they should sub- 
sequently accrete. 

At a higher level of sophistication, sink properties 
(masses, positions, velocities, accretion rates) are used to in- 
fer stellar mass functions, binary statistics, cluster kinemat- 
ics and the star formation rate in turbulent molecular clouds 

20021: 



(e.g. iBonnell et al l l200ll: iKitsionas fc Whitworth 



[ Bate et all I2QQ3I: iGoodwin et al.1 l2004allbl: IBonnell et al l 
l2004l:lGoodwin et a l. 2006; Bonnell & Ba te 2006; Bate 20091: 
[ Attwood et al.1 120091: IStamatellos et al.1 l201ll : IWalch et all 
12012b iFederrath fe Klessenl l2012h . Thev can also be 
used to model radiative and mechanical feedback from 
protostars, using simple models or phenomenological 
prescriptions for the physics inside the sink (e.g. 
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Stamatel los et all [20051: bale et all 120051; iKrumholz etaD 



2007: IStamatellos fe Whitworth 2009b; Offner et al.1 Ig009: 

et al.1 

20101: iDale fc Bonnelll 120111; IPeters et al l l201ll : lBatdl201S 
Cunningham et al. 201*3: bale Bo nnell 201§). 



In cosmological simulations, sink particles are used 
to model the format ion and growth of black holes, and 
hence AGN f eedback ([Springel et alJl2005l: iDi Matteo etafl 



2011 



2005, 120081: Johanss on et all I2OO9I : IPebuhr et all [2 010, 



Choi et al 



l2Q12h . For example, in Springel et alJ 
(|2OO50 the efficiency of AGN fee dback is estimated using 
the Bondi-Hoyle-Littleton theory (|Hovle &; Lvttleton| ll939: 
iBondi fe Ho^y|l944l : lBondilll952h acting on scales of a few 
tens of parsecs, whereas the actual Schwarzschild radius of a 
typical supermassive black hole is of order a few solar radii. 

Two main algorithms are currently used to mode l star 
form a tion, Smoothed Particle Hy drodynamics (SPH; iLucvl 
Il977l : iGingold &; Monaghanlll977h and A d aptive Mesh Re- 
finement (AMR: IBerger fe Colellal Il989l: IBerger &; Oliger] 



1 19841 : iDezeeuw & Powelll ll993: MacNeice et al [[2000|). This 
paper is concerned with sinks in SPH. 

In SPH simulations of star formation, two require- 
ments are critical. First, the mass of a single SPH 
parti cle, 



s hould be chosen (e.g. [B ate & Bu rkertl 



19971: IWhitworthl Il998l: lAttwood et al.1 l2009l : iBatd 120091 : 

Stamatellos &; Whitworth] 2009a ). or adjusted (e.g. 



Kitsionas &; Whitworth! 120021 ). so that the Jeans mass 



is always resolved. Second, sink creation should be consid- 
ered only if the density exceeds a user-defined threshold, 
Psink 5 an d Psink should in turn exceed the density at which 
the minimum Jeans mass obtains, so that the code can 
capture properly the full range of fragmentation, and the 
formation of first cores. Thus our main concern here is 
with sinks created from gas that is already being heated by 
adiabatic compression; for a discussion of how to treat sinks 
created at lower densities fro m approximatel y isoth ermal 
gas the reader is refered to iFederrath et al.l (|2010h . For 
contemporary local star formation, involving molecular gas 
at T ~ 10 K, these requirements reduce to m SPH < 10 -5 M 
and p SINK ^> 10 -13 gem -3 . Once the density exceeds p SINK , 
the criteria that are applied to determine (i) whether a 
sink actually is created, (ii) what matter it immediately 
assimilates, and (iii) what matter it subsequently accretes, 
must be formulated so as to obtain close correspondence 
between the properties of the sinks that do form, and the 
stars that should have formed, under the given physical 
circumstances. 

In AMR simulations of star formation, the computa- 
tional domain is divided adaptively into finer meshes - and 
as appropriate these meshes are subsequently de-constructed 
- so as to deliver high resolution only when and where it 
is needed. There is a minimum mesh spacing, ^ MIN , below 
which refinement is prohibited. I n the original implem en- 
tation of sink particles in AMR ([Krumholz et al.ll2004l ). a 
sink is introduced wherever the local Jeans length falls be- 
low 4^ MIN and is therefore no longer adequately resolved. 
There are then rules to determine the rate at which mass 
is accreted onto the sink from the surrounding cells, and to 
identify the circumstances under which neighbouring sinks 
are merged. The angular momentum of the acc reted mass is 
left behind in the cell from which it is accreted. Offner et al. 
(2008) have shown that this implementation is quite ro- 



bust, in the sense that there is good correspondence be- 
tween the properties of sinks formed in simulations having 
the same initial conditions but different £ MIN . To capture the 
full range of fragmentation, and the formation of first cores, 
requires ^ MIN <C 3AU, for contemporary local star forma- 
tion. We note that in some AMR sink implementations (e.g. 
IFederrath et al.ll20ld ) the angular momentum of matter ac- 
creted by the sink is not left behind, but is assimilated by 

the sink. 

Recently, IFederrath et al ] (|2Q10h have simulated the for- 
mation of a star cluster with SPH and AMR, using compa- 
rable resolution, and have developed a set of sink-creation 
criteria which ensure that the two methods give very simi- 
lar results, in terms of the masses, locations and accretion 
histories of the sinks formed. However, although these simu- 
lations investigate the consequences of different sink creation 
criteria, they do not explore in detail the dependence of the 
sink properties on the resolution, or on the user-defined pa- 
rameters of sink creation. The present paper addresses these 
issues. 

In Section [2] we describe in detail a new algorithm for 
the creation and evolution of sinks, explaining the motiva- 
tion for each of its features; we call these NewSinks. In 
Sections [3] and [4] we describe the sink algorithms that have 
been used previously in SPH simulations of star formation - 
hereafter standard sinks - focussing on the particular imple- 
mentations that we use for comparison tests. Since all but 
one of the new features in NewSinks relate to their evolu- 
tion, rather than their creation, Section [3] defines OldSinks, 
which are created in the same way as NewSinks, but evolved 
in the manner of standard sinks, thereby allowing us to 
isolate those aspects of NewSinks that have to do purely 
with sink evolution. On the other hand, since most extant 
SPH simulations of star formation that invoke sinks actually 
use a less robust algorithm to create sinks than that used 
by NewSinks and OldSinks, Section g] defines UrSinks, 
which are created and evolved in the manner of standard 
sinks. The basic properties of the different types of sink are 
summarised in Table [1] Sect ion [5] presents the results of tests 
performed with the different types of sink to demonstrate 
that when NewSinks are used the results obtained are only 
weakly dependent on the resolution and the user-defined pa- 
rameters of sink-creation, whereas when OldSinks are used 
the results are strongly dependent and not converged, and 
when UrSinks are used the results are divergent. General 
aspects of the different types of sink are discussed in Section 
[6] Our main conclusions are summarised in Section 



2 THE NewSink ALGORITHM 

Here we describe and justify the procedures used to define 
the internal structure and translation of a NewSink (Section 
12. ip . to trigger the creation of a NewSink (Section [272]) , to 
accrete matter onto an existing NewSink (Section 12. 3[) , to 
update the properties of a NewSink (Section 12. 4[) , and to 
pass the angular momentum acquired by a NewSink back 
to the surrounding matter (Section |2.5|) . 
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Type of Sink 




NewSink 


OldSink 


UrSink 


Creation criteria: 


Equation: 










density 


Pi > Psink 




* 


* 


* 


overlap 


|r; - r a | >X SINK hi + R s 


(3) 


* 


* 


* 


potential minimum 


(pi < MIN {</>j } 


CO 


* 


* 




Hill Sphere 


P^>Phill= ••• 


© 


* 


* 




acceleration divergence 


(V-a)i<0 


([33]) 






* 


velocity divergence 


(V-v)i<0 


(f34| 






* 


non-thermal energy 


£ NT <0 


(1351) 






* 


Structure: 

interaction-zone 
exclusion-zone 


Evolution: 

regulated accretion 
instantaneous accretion 




■ 


■ 


■ 



Table 1. The creation criteria, intrinsic structures, and evolution procedures for the different types of sink. Column 1 gives the various 
creation criteria and the two different intrinsic structures. Column 2 gives equations for the various creation criteria, and the two different 
modes of evolution. Column 3 gives equation numbers. In columns 4 (NewSink), 5 (OldSink) and 6 (UrSink), stars indicate which 
creation criteria are invoked, and filled squares indicate which structure/evolution combination is used. Relative to NewSinks, OldSinks 
have the same creation criteria, but different evolution procedures; UrSinks have different creation criteria and different evolution 
procedures. Most extant SPH simulations of star formation have been performed with standard sinks like UrSinks. 



2.1 Internal structure and advection of a NewSink 



Pi > PsiNK ' 



(2) 



A NewSink comprises a central point mass (hereafter, the 
point-mass) and a concentric spherical volume of radius R s 
(hereafter, the interaction- zone) that moves with, and re- 
mains centred on, the point-mass. When we refer to the 
mass, M s , position, r s , velocity, v s , or angular momentum, 
L s , of a NewSink, we mean the mass, position, velocity and 
spin angular momentum of its point-mass. Once created, a 
NewSink is tracked with the same integration routine as an 
SPH particle, modulo that it only experiences gravitational 
acceleration. The interaction-zone contains live SPH parti- 
cles that have not yet been assimilated by the point-mass, 
and which are intended to represent the continuation of any 
accretion disc, inflow stream, or other coherent structure im- 
mediately outside the interaction-zone. The interaction-zone 
is one of the key features that distinguishes a NewSink from 
a standard sink. 



2.2 Criteria for creation of a NewSink 

The creation of a NewSink is triggered by an SPH parti- 
cle. Specifically, if SPH particle z, having smoothing length 
hi, satisfies the four criteria below, it is replaced with a 
NewSink, s. At the moment of creation, the point-mass of 
s has the same mass, position and velocity as i did, and the 
radius of its interaction-zone is 



R s — X slNK hi 



(i) 



X SINK is a user defined parameter, chosen so that the neigh- 
bours of i all fall inside the NewSink at the moment of its 
creation. Thus, for the M4 kernel, we advocate ^ SINK = 2. 
By adopting a larger X SINK one obtains smoother accretion 
onto a NewSink, but at the expense of coarser resolution. 



2.2.1 Density criterion 

The first criterion for SPH particle i to trigger the creation of 
a NewSink is that the density at its location, pi = p SPH (r^), 
should exceed a user-defined threshold, p SINK , i.e. 



We advocate p SINK ^ 1CT 11 gem -3 , so that NewSinks nor- 
mally form in condensations that are well into their Kelvin- 
Helmholtz contraction phase. 



2.2.2 Sink-overlap criterion 

The second criterion for SPH particle i to trigger the forma- 
tion of a NewSink is that it should not, at the moment of 
creation, overlap another NewSink, i.e. 

|r< - r a / 1 > X slNK hi + R s > , (3) 

for all pre-existing NewSinks, s' . 

2.2.3 Gravitational potential minimum criterion 

The third criterion for SPH particle i to trigger the forma- 
tion of a NewSink is that its gravitational potential, 0i, 
should be lower than that of all it neighbours, j, i.e. 

4>i < MIN {<f)j} . (4) 

This criterion, introduced by iFederrath et al ] (|201Ch . en- 
sures that NewSinks are only created near resolved po- 
tential minima (and hence, presumably, resolved density 
peaks). 

2.2.4 Hill Sphere criterion 

The fourth criterion for SPH particle i to trigger the forma- 
tion of a NewSink is that its density should satisfy 



Pi > Phi 



3X HILL (-Ar^-Aa^) 
4^G|Ar^| 2 



(5) 



for all pre-existing NewSinks, s'. Here AT HILL is a user- 
defined parameter with default value AT HILL = 4, Ar is f = 
Ti—r s / and Aa is / = a — a s / are the position and acceleration 
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of particle i relative to s' This criterion deals with the sit- 
uation where a condensation undergoing Kelvin-Helmholtz 
contraction has a NewSink s' in its interior, but is much 
more extended than that NewSink. Eqn. ([5]) ensures that a 
second NewSink, seeded by SPH particle i, can form in the 
outskirts of the condensation only if there is a density peak 
at Yi that dominates the local gravitational field. 

2.3 Criteria for accretion by a NewSink 

2.3.1 The SPH particles inside a NewSink 

Under normal circumstances, SPH particles that enter the 
interaction-zone of a NewSink are not accreted immedi- 
ately by its point-mass. In the first instance, they are simply 
added to a list of the SPH particles inside the interaction- 
zone, its interaction- list. These SPH particles may subse- 
quently be accreted by the point-mass, possibly over sev- 
eral timesteps, but they may leave the interaction-zone be- 
fore this happens. If an SPH particle finds itself inside the 
interaction-zones of more than one existing NewSink, it is 
added to the interaction-list of the NewSink whose point- 
mass is closest. In this context, it is appropriate to note that, 
although a NewSink is never created overlapping another 
NewSink, its motion may there after lead to an overlap . We 
do not allow sinks to merge (cf. iKrumholz etaD 12004 ). 

To determine whether, and how quickly, the SPH par- 
ticles in the interaction-zone of an existing NewSink are 
assimilated by its point-mass, we consider the limiting cases 
of spherically-symmetric radial accretion, and disc accretion. 



2.3.2 Timescale for spherically symmetric radial accretion 

In spherical symmetry, the rate at which mass flows inwards 
across a spherical surface at radius r is 



M(r 



4vrr p(r)v RAU (r) . 



(6) 



At the position of an SPH particle j, in the interaction-zone 
of NewSink s, this inflow rate becomes 



47r | Ar js | Ar js • Av js p 3 - . 



(7) 



The timescale for radial accretion is obtained by dividing 
the net mass of the SPH particles in the interaction-zone by 
a weighted sum of the inflow rates: 



<*rad > s 4n ^{ | Ai> | Ar js • Av JJ m J W( | Ar js \ , H s ) } ' 



W = £{m J W(|Ar J .|,.ff.)M} 



(8) 
(9) 



The weighting here uses the kernel- function, W. The 
smoothing length of the sink, H s , is adjusted so that the ex- 
tent of the kernel function equals R s ; for example, with the 

1 Here, and in the sequel, we adopt the convention that Aq a ^ = 
Qa— qb is the difference in some quantity q between two particles a 
and b, evaluated at the same time; thus, for example, Ar^ = — 
Yj is the instantaneous position of SPH particle i relative to SPH 
particle j. In contrast, we adopt the convention that 5q c is the 
change in some quantity q c associated with particle c, due to the 
evolution of the system; thus, for example, (5L S is the increment 
to the angular momentum of sink s, and 5t s is the timestep of s. 



M4 kernel, H s = R s /2. Consequently SPH particles in the 
outer regions of the interaction-zone make a smaller contri- 
bution than those closer to the point-mass. W ensures that 
the sum is accurately normalised. 

2. 3. 3 Timescale for disc accretion 

We m odel disc accretion using the IShakura &; Sunvaevl 
(| 19731 ) prescription, which conflates all possible angular 
momentum transport mechanisms into a single parameter, 
a ss . For a low-mass disc in approximate Keplerian rotation 
around a star of mass , the accretion timescale at radius 
R is then ~ (GM+R) 1 / 2 a~ 2 , where a is the local sound 
speed. We therefore compute a kernel- weighted mean over 
all the SPH particles in the interaction-zone, 

(GM S ^ 1/2 ■ 



a ss W 



-E 



Ar Js ^ 2 



m 3 W(\AY 3 



Pjd 2 



(.10) 



Both observa tional estimates of protostellar accretion disc 
lifeti mes (e.g. Hartmann 19981). and theoretical simulations 
(e.g. iForgan et al.l 201(? V. suggest that a ss = 0.01 is a suit- 
able default value, but there are large uncertainties. 



2.3.4 Net timescale for accretion 

The timescale for accretion onto the point-mass is given by 

(id 



t = (t (+ ) f 

L ACC VRAD/s VDISC/s 



/ 



MIN • 



\E n 



(12) 



Here ^ ROT an d ^ GRAV are - respectively - the net rota- 
tional and gravitational energies of the SPH particles in the 
interaction-zone, relative to the point-mass, i.e. 



(13) 



2£,-{m i |Ai > .L INT P}' 
GM S 



^{mjAr Js xAv JS } 



?>)■ 

(14) 
(15) 



is the net angular momentum of the SPH particles in the 
interaction-zone, relative to the point-mass, and is a func- 
tion that encapsulates the kern el-smoothing of the gravita- 
tional potential (see Eqn. 15 in lHubber et aljfeoilh . 

Eqn. [11] is adopted because it gives the correct limiting 
behaviour. If the SPH particles inside the sink are in rota- 
tional equilibrium, t ACC — > (t DISC ) s . Conversely, if they are 
not rotating at all, t ACC — > (t RAD ) s . 

2. 3. 5 Excising SPH particles 

Normally, the mass accreted by the point-mass at the end 
of the current timestep, (£,£ + St s ), is 



5M ACC = M INT 



exp 



St 3 



(16) 
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In the first instance, this mass is removed from the SPH 
particle closest to the point-mass. If the mass of this SPH 
particle is less than 5M ACC , the remainder is removed from 
the second closest SPH particle, and so on, until the whole of 
SM ACC has been removed; because this only affects particles 
close to the point-mass, it is extremely rare that particles 
with reduced mass leave the interaction-zone. There are two 
circumstances under which the procedure outlined above is 
superseded. 

(i) If the total mass of SPH particles in the interaction- 
zone of a NewSink presently exceeds the total mass of SPH 
particles inside the interaction-zone at the time it was cre- 
ated, M MAX , then t ACC is decreased artificially by a factor 
t ACC —> t ACC /(M INT /M MAX ) 2 , in order to accrete the excess 
mass more rapidly. This is similar to the procedure used to 
transfer mass from nea rby grid cells to sink particles in AMR 
(iKrumholz et alJl2004h . 

(ii) If the timestep, Stj, for an SPH particle, j, 
in the interaction-zone of NewSink s, satisfies Stj < 
7 S (Rg/GMs) 1 / 2 , j is immediately accreted - in its total- 
ity - by s. 7 S is a tolerance parameter, with default value 
0.01. This device is intended to moderate the situation where 
the SPH particles inside the interaction-zone have formed a 
Toomre unstable disc. The assumption is that this will lead 
to efficient transport of angular momentum by gravitational 
torques, and hence very rapid inspiral onto the point-mass. 

This piecemeal leaching of mass, from the SPH particles 
in the interaction-zone, onto the point-mass at its centre, is 
termed regulated accretion, and is a critical feature of the 
NewSink algorithm. It ensures that, as long as the sink is 
accreting, there are SPH particles in the interaction-zone, 
and therefore there are no very steep gradients across the 
boundary of the sink. 

The idea of transferring mass piecemeal from an SPH 
particle to a sink is not new. T his device has been im- 
plemented by lAnzer et al ] <|l987h , in simulating the accre- 
tion of a supersonic wind by a neutron star, and has sub- 
sequently been used i n the same context by, for example, 
iBoffin fe Anzerl (|l994h . In their algorithm, mass is trans- 
ferred to a pre-existing sink, from all the particles in the 
computational domain - but at a rate that is strongly 
weighted towards those presently near the sink; the weight- 
ing function has to be tuned to produce an acceptable rate 
of accretion onto the sink. However, here we are consider- 
ing dynamically created sinks and subsonic ambient flows. 
Moreover, in our algorithm, the particles that leach mass 
to the point-mass are always the ones near the centre of 
the interaction-zone, which never leave the interaction-zone; 
consequently, we only have to contend with a small number 
of particles having time- varying mass, and none outside of 
the interaction-zone. 



r' s = M'.- 1 (M.r. + '52{6m j T j }J , 

L' s = L s + M s Ar ss , x Av ss , 

+ ^2 {5m j Ar js , x Av jV } • 



(18) 
(19) 

(20) 



Here, the summation is over all SPH particles, j, that lose 
mass, 5m j, to the point-mass, and the updated values are 
denoted by primes. Once an SPH particle has zero mass, it 
is removed from the simulation altogether. 



2.5 Angular momentum feedback from a NewSink 

If the point-mass were to assimilate all the angular momen- 
tum of the SPH particles it accreted, as implicit in Eqn.[20l 
it would spin so fast that it could not reach stellar densi- 
ties. This is what happens with all standard sinks in SPH: 
they are sinks of both mass and angular momentum, and 
this is highly unrealistic. In reality, if the material inflowing 
towards a protostar has high specific angular momentum, it 
first falls onto an accretion disc, and then spirals inwards by 
transferring the bulk of its angular momentum to other ma- 
terial further out in the disc or envelope. Therefore, at the 
end of each timestep, we reduce the angular momentum of 
the point-mass by transferring some of its angular momen- 
tum to the SPH particles in the surrounding interaction- 
zone. Since this transfer of angular momentum is presumed 
to be effected by viscous torques in an accretion disc, the 
amount of angular momentum transferred is given by 



\SL S 



1 — exp 



St s 



(21) 



(^DISc) / 

We achieve this by giving each SPH particle, j, in the 
interaction-zone of s an impulse of velocity 

\SL S \ L s xAr js 



V • {mj Ar JS xLs x Ar^s} 



(22) 



Thus each SPH particle in the interaction-zone of s receives 
an impulse of velocity proportional to its distance from the 
rotation axis of the point-mass. 

In order to compensate for these impulses of velocity, 
the point-mass receives impulses of momentum and angular 
momentum given by 



3 

5~L S — — {rrij Arj S x5vj} . 



(23) 
(24) 



2.4 Updating the properties of a NewSink 

At the end of each timestep, St s , the mass, M s , position, r s , 
velocity, v s , and angular momentum, L s , of the point-mass 
are updated to take account of accretion, 



(17) 



As a consequence, the net angular momentum (invested in 
SPH particle motions and the spins of point-masses) is accu- 
rately conserved. At the same time, the amount of angular 
momentum invested in the spins of point-masses is small, 
both because it is continually transferred to the SPH par- 
ticles in the interaction-zone, and because the close-in SPH 
particles from which it assimilates mass have normally had 
to lose a lot of angular momentum to get close-in in the first 
place. This is a critical feature of the NewSink algorithm. 
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It ensures that NewSinks do not act as sinks for angular 
momentum, and that the point-mass can only grow in mass 
at a rate influenced by how fast the close-in SPH particles 
are able to transfer their angular momentum to other SPH 
particles further out in the interaction-zone and beyond. 



3 THE OldSink ALGORITHM 

Many variants of standard sink have bee n used previously 
in SPH simulations of star formation (e.g. [ Bate et al.l ll995; 
iFederrath et al.1 l201(]l : IWadslev etHI l201lh . and~it is not 
practical to present test results for all of them. Here, and 
in Section 2J we describe the particular representative im- 
plementations that we use for comparison tests. The choices 
have been made with a view to distinguishing the problems 
associated with sink creation from those associated with sink 
evolution. In this Section we define OldSinks, which are cre- 
ated in the same way as NewSinks, and therefore avoid the 
problems normally associated with the creation of standard 
sinks. This enables us to isolate the problems associated with 
the evolution of standard sinks, viz. that they all appear to 
seriously corrupt the hydrodynamics near the sink bound- 
ary (by introducing excessively steep gradients), to act as 
sinks for angular momentum, and to produce results that 
are strongly dependent on the user-defined parameters of 
sink creation and evolution. 



3.1 Internal structure and advection of an 

OldSink 

An OldSink comprises a central point mass (the point- 
mass) and a comoving concentric spherical volume of radius 
R s (hereafter the exclusion- zone). The mass, M s , position, 
r s , velocity, v s , and angular momentum, L s , of an OldSink 
refer to its point-mass. The exclusion-zone, as its name im- 
plies, generally contains very few SPH particles - and this 
is one of the critical features that distinguishes an OldSink 
from a NewSink. Once created, an OldSink is tracked with 
the same integration routine as an SPH particle. 



3.2 Criteria for creation of an OldSink 

The criteria for creation of an OldSink are the same as 
for creation of a NewSink, viz. that an SPH particle, i, 
has density exceeding p SINK (Eqn. [2}, that the resulting 
sink does not - at its creation - overlap a pre-existing sink 
(Eqn. [3}, that i has lower gravitational potential than all 
its neighbours (Eqn. |4}, and that i satisfy the Hill Sphere 
criterion (Eqn. [5}. If an SPH particle i satisfies these four cri- 
teria, an OldSink is created with the same mass, M s = rai, 
position, r s = n, velocity, v s = v*, and angular momen- 
tum, L s = 0, as i. The radius of the exclusion-zone is set 
to R s — X SINK hi, where ^ SINK is a user-defined param- 
eter. With the M4 smoothing kernel, the default value is 
X SINK = 2, so that the exclusion-zone is the same as the 
smoothing volume of i. 



3.3 Criteria for accretion by an OldSink 

An SPH particle j is assimilated by the point-mass of an 
OldSink, s, if it satisfies two criteria. 



First, it must have entered the exclusion-zone of s. i.e. 



< R s 



(25) 



Second, the mutual non-thermal energy of j and s (i.e. 
their mutual bulk-kinetic plus gravitational energy), 

|2 



2(mj + Ms) 



■ + - 



GrrijMs 



Ar, 



hi 



Ar, 



(26) 



must be negative. If more than one OldSink is minded to 
accrete j, j is accreted by the one whose point-mass is clos- 
est. When an SPH particle i triggers the formation of an 
OldSink, s, these conditions are normally satisfied by all, 
or most, of the neighbours of i, so s immediately assimilates 
all, or most, of t he ne ighbours of i. 

iBate et all |l995) use two additional criteria to deter- 
mine whether an SPH particle, j is assimilated by a stan- 
dard sink, s. First, the angular momentum of j relative to 
the point-mass of s should be less than it would be if j were 
in a circular orbit of radius R s about the point-mass, i.e. 



< (GM 8 R t 



\V2 



(27) 



Second, if an SPH particle j passes within a small distance 
(O.IRs or less) of the point-mass, all other criteria are ig- 
nored, and it is assimilated by the point-mass. We find that 
thes e criteria d o not h ave a significant effect on the results. 

IBate et all (|l995h have also developed a procedure for 
extrapolating the density and velocity fields outside a sink, 
and using this information to generate correction terms to 
account for the missing SPH particles inside the exclusion- 
zone. However, these terms ha ve not been includ ed in sub- 
sequent implementations of the lBate et al.l (|l995h algorithm 
(Bate, private communication), and we also do not include 
them in OldSinks (or UrSinks, see Section [4}. 



3.4 Updating the properties of an OldSink 

Accretion of an SPH particle by an OldSink, s, is imme- 
diate and complete; the SPH particle is removed from the 
simulation, and the properties of s are updated according to 



M' s = Ms + Y^imj}, 

3 

vi = M' s _1 ^M s v s + ^ {rrij Vj}^ , 

L' s = L 3 + M s Ar ss ,xAv ss . 

+ ^2{ m j AryxAvy} 



(28) 
(29) 
(30) 

(31) 



the primed variables represent the new properties of s, and 
the sums are over all the SPH particles j accreted during 
that timestep. We term this instantaneous accretion. 

We note that, if an OldSink grows by disc accretion, 
most of the SPH particles that it accretes only just satisfy 
Eqn. 1251 and so its specific angular momentum is 



|L,| in 20 2 -1 f M s 

^—r 1 ~ 10 cm s — — 

Ms \M 



1/2 
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( m SPH A / f PsiNK \ ^ /00\ 

Vio- 5 M y v i0_13 g cm " 3 / 

This is so large that the point-mass can not condense 
to stellar densities (for example, the Sun spinning at 
break-up speed only has specific angular momentum ~ 
cm s ), unless m srH is very small 
(impractically high mass-resolution) and/or p SINK is very 
high (so high as to negate the advantages of using sinks). 



4-2.3 Vv criterion 

The fourth criterion for SPH particle i to trigger the creation 
of an UrSink is that the divergence of the velocity at the 
position of i should be negative, 

(V-v), < 0, (34) 

otherwise the density there must be decreasing. This crite- 
rion is used by Wadslev et al ] (1201 ih . 



4 THE UrSink ALGORITHM 

In this Section we define UrSinks, which are created and 
evolved in the manner of standard sinks, and are represen- 
tative of the sinks that have been used in almost all extant 
SPH simulations of star formation. 



4.1 Internal structure and advection of an UrSink 

An UrSink has the same internal structure as an OldSink, 
viz. a central point-mass surrounded by a comoving concen- 
tric spherical exclusion-zone of radius R s . The point-mass 
has mass, M s , position, r s , velocity, v s and angular mo- 
mentum L s , and is advected like an SPH particle. 



4.2 Criteria for creation of an UrSink 

Some of the criteria for creation of an UrSink are different 
from those used for NewSinks and OldSinks. If an SPH 
particle i satisfies the five criteria detailed below, an UrSink 
is created with the same mass, M s = im, position, r s = r*, 
velocity, v s = v^, and angular momentum, L s = 0, as i. 
The radius of the exclusion-zone is set to R 3 — X slNK hi, 
and if the M4 smoothing kernel is used, the default value 
is A SINK — 2, so that the exclusion-zone is the same as the 
smoothing volume of i. 



4-2.1 Density and overlap criteria 

The first two criteria for SPH particle i to trigger the cre- 
ation of an UrSink are the same as for NewSinks and 
OldSinks, viz. that its density should exceed a user-defined 
threshold, and that the resulting UrSink should not overlap 
a pre-existing UrSink (see Sections [2.2.11 and f2.2.2[) . 



4.2.2 Va criterion 

The third criterion for SPH particle i to trigger the creation 
of an UrSink is that the divergence of the acceleration at 
the position of i should be negative, 

(V-a), < 0, (33) 

otherwise the possibility exists that the dense g as is about 
to be torn apart tidally. T h is crit erion is used bv lBate et al.l 
(|l995h and lWadslev et all (|201lh . 



4.2.4 Non-thermal energy criterion 

The fifth and final criterion for SPH particle i to trigger the 
creation of an UrSink is that the net non-thermal energy of 
i and its neighbours j, in the centre-of-mass frame, should 
be negative, 

E NT = E GRAY + £ KIN < . (35) 
Here, 

(36) 

is the self-gravitational potential energy, 

E Km = ^{m.lvi-v/} (37) 

3 

is the net kinetic energy, and Vj = ^2 {rrijVj} /^2{rrij} is 
the centre-of-mass velocity. Here all the sums over the neigh- 
bours j include i itself. This criterion is somewhat weaker 
than the energy criteria used by other standard algorithms 
(see Section 14.2.51 below) . 

4.2.5 Alternative creation criteria 

In this Section, we list - for completeness - some of the 
additional criteria used in other standard sink creation al- 
gori thms, but not for UrSi nks. 

bate et all (|l995h and IWadslev et~aT1 (|201lh apply ad- 
ditional energy criteria for creating a standard sink, viz. 



^THERM + 2 ^GRAV < ^ , (38) 

^ROT ~^ ^THERM ~^ ^GRAV ^ ^ J (^9) 

-^KIN + ^THERM + ^GRAV < ^ 1 (^0) 

Here 

^THERM | ( 7j _ 1) } ( 41 ) 



is the net thermal energy (aj and jj are the isothermal sound 
speed and ratio of specific heats for SPH particle j), 

J? I^INT I (ac)\ 

2£,{m,|Ar, 3 .L INT P} ^ 
is the net rotational energy, 

Lint {™j ( r J ~ *3 ) X ( V J ~ Vj ) } (43) 

3 

is the net angular momentum, and fj = { m j r j} I { m i} 
is the centre of mass. All the sums over the neighbours j in- 
clude i itself. Although these energy criteria are somewhat 
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stiffer than Eqn. ([35)) . they do not produce significantly dif- 

fere nt reaults. 

iFederrath et all (|2010h and IWadslev etHI (|201lh also 
use the potential-minimum criterion (i.e. Eqn. [4j), and in 
addition they advocate checking that the flow is convergent 
in all directions, by computing the eigenvalues of dvi/dxj 
(here i and j are identifiers for the Cartesian coordinates) 
and requiring them to be individually negative. However, 
these criteria have not been implemented in the majority of 
extant SPH simulations of star formation. 



4.3 Criteria for accretion by an UrSink 

An SPH particle j is immediately accreted by the point-mass 
of an UrSink s, if it falls within the exclusion-zone of s, and 
their mutual non-thermal energy is negative, just as for an 
OldSink fEqns. [25l I26[) . Under normal circumstances, a 
newly-created UrSink immediately assimilates most of the 
neighbours of the SPH particle that triggered its creation. 



4.4 Updating the properties of an UrSink 

Each time an UrSink assimilates one or more SPH parti- 
cles, the properties of the UrSink are updated according to 
Eqns. ([28]) through ([3T]) . just as for an OldSink. As with 
an OldSink, this often leads to the unphysical acquisition 
of large amounts of angular momentum. 



5 TESTS 

We have used the following four tests to compare NewSinks 
with OldSinks and UrSinks, and to explore the effect of 
changing the user-prescribed parameters for sink creation 
(Psink an d -^sink) an d the SPH resolution (A/" SPH ): isother- 
mal Bondi accretion (Section 15. 1[) . the collapse of a rotating 
Bonnor-Ebert sphere (|5.2[) . the Boss-Bodenheimer test ([5.3[) , 
and the fragmentation of a turbulent prestellar core (|5.4[) . 
Parameter values used for the different tests are summarised 
in Table 2. Although all four tests have been performed with 
all three types of sink, we do not present and discuss results 
obtained with UrSinks for the first three tests; this is be- 
cause the results obtained with UrSinks are very chaotic 
and show no sign of converging; the turbulent prestellar core 
test suffices to demonstrate this. 

All the tests are pe rformed with the SE REN SPH code, 
using its default options (|Hubber et al.ll201lh . We set 77= 1.2, 
so that the smoothing-length of SPH particle i is hi = 
1.2(ra SPH / pi) 1 ^ 3 and the mean number of neighbours is 



a/; 



^58. 



5.1 Bondi accretion 

The first test models Bondi accreti on. The accre tion rate for 
the transsonic isothermal solution (jBondilll952h 



e 3/2 7T G 2 Ml p c 



(44) 



where e is the base of natural logarithms, M* is the mass 
of the central star, p Q is the density at infinity, and a Q 
is the isothermal sound speed. The solution is obtained by 
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Figure 1. The normalised accretion rate (M s / M-qo) as a func- 
tion of the normalised sink radius (R s /R SONIC ) for isothermal 
Bondi accretion. Results obtained with NewSinks are shown as 
red diamonds connected by a solid line, and those obtained with 
OldSinks are shown as black stars connected by a dashed line. 



neglecting the self- gravity of the inflowing gas. A critical role 
is played by the sonic radius at 

GM* 



R 



2al 



(45) 



For r ^> R SONIC ^ ne inward flow is subsonic and the pressure 
acceleration (— p~ x VP) plays an important role. Conversely, 
for r <C ^sonic the inward flow approaches free fall and the 
pressure acceleration is unimportant. 

The density and velocity profiles for the transsonic solu- 
tion are obtained by solving the Bernouilli equation numeri- 
cally. The initial conditions are set up by relaxing a periodic 
cube to produce a uniform-density glass, cutting a sphere 
containing 5 x 10 5 particles from this cube, stretching the 
sphere to produce the required density field, and giving each 
SPH particle the inward radial velocity appropriate to its po- 
sition. The outer boundary of the sphere is at ~ 20 R SON i C , 
and the simulations are followed to time t END = 2GM*/a 3 J . 

This test does not involve sink creation, and therefore 
we do not specify p SINK or X SINK . Instead, a sink is placed at 
the origin from the outset. Its point-mass has mass M s = M 
and, although it accretes SPH particles and has an accretion 
rate, M s , its mass is held fixed. This is self- consistent, since, 
if the accreting gas is sufficiently rarefied to have negligible 
self-gravity (compared with the gravity of the central star), 
the increase in the central mass by t END is negligible. The 
purpose of this test is to evaluate the treatment of radial 
accretion. 

The system is evolved with a NewSink or an OldSink, 
with the sink radius set to different multiples of the sonic 
radius, R s /R SONlc = 1/8, 1/4, 1/2, 1, 2, 4, 8. By t END , the 
accretion rate is steady, and the rarefaction wave propagat- 
ing in from the outer boundary is still in the outer parts 
of the computational domain. Fig. [1] shows how the accre- 
tion rates depend on ^s/^ SO nic- Since this test does not 
involve sink creation, the results obtained with an UrSink 
are exactly the same as those obtained with an OldSink. 

As long as R s < R SONIC , both the NewSink and the 
OldSink model Bondi accretion well, giving accretion rates 
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Test — ► 


Default 


Bondi Accretion 


Rotating Bonnor-Ebert 


Boss-Bodenheimer 


Turbulent Core 


PsiN K /(g cm " 3 ) 


10- 11 




io- 11 


io- 13 , io- 12 , io- 11 


io- 11 , io- 10 , io- 9 


^SINK 


2 




2 


2, 4 


2 




4 




4 


4 


4 


a ss 


0.01 


0.01 


0.01 


0.01 


0.01 


Rs / R^ONIC 




iiil, 2, 4,8 








NewSink 




* 


* 


* 


* 


OldSink 




* 


* 


* 


* 


UrSink 




(*) 


(*) 


(*) 


* 



Table 2. Default parameter values and values used for the different tests. Row 1 gives the name of the test. Rows 2 to 6 give the 
parameter values, and rows 7 to 9 indicate which types of sink were tested, with brackets indicating the tests that are not discussed. 



~ 10% below the analytic value. Under this circumstance, 
the material flowing into the sink is close to free fall and 
so the pressure acceleration is unimportant; the fact that 
OldSinks give rise to inaccurate evaluations of the pressure 
and viscous forces near R s is then of little consequence. 

However, for R s > R SONlc , an OldSink seriously over- 
estimates the accretion rate. Under this circumstance, the 
pressure acceleration is an important factor controlling the 
flow of material into the sink, and the steep outward pres- 
sure gradient at the edge of an OldSink artificially increases 
the accretion rate; for Rs/R SONIC = 8, M s is ten times the 
analytic value. In contrast, with a NewSink there is a small 
increase in the accretion rate with increasing R s , but, even 
for R s /R SONlc = 8, M s is within ~ 1% of the analytic value. 

Although this test only deals with one possible accre- 
tion mode, all spherical accretion modes are likely to suffer 
from the same problem, if the inflow is subsonic near the sink 
boundary. This is a serious concern in simulations of star for- 
mation that set p SINK high, in order to follow protostars well 
into their Kelvin-Helmholtz contraction phase. The flow of 
gas into the sink is then subsonic, and an OldSink will ex- 
perience artificially enhanced accretion. A NewSink should 
capture the accretion rate much more accurately. 

5.2 Rotating Bonnor-Ebert sphere 

The main mode of accretion in star formation involves discs. 
The accretion rate is then controlled by the torques that re- 
distribute angular momentum. There are no analytic solu- 
tions for this mode, so testing is less straightforward. The 
test we perform involves a rigidly rotating spherical gas 
cloud having mass M , initial boundary radius 10 3 AU, and 
initial angular speed 3.54x 10 -12 s _1 . The gas is and remains 
isothermal with sound speed 0.20kms _1 (corresponding to 
molecular gas at T= 11 K). However, the main purpose of 
this test is not to explore sink creation from isothermal gas, 
but rather to evaluate the treatment of disc accretion. 

The initial density profile is that of a critical Bonnor- 
Ebert Sphere, i.e. one with dimensionless boundary radius 
£ B = 6.45, but the cloud is not in equilibrium. The initial 
ratios of thermal and rotational energy to gravitational en- 
ergy are a = 0.09 and f3 = 0.59. The cloud is modelled with 
10 5 SPH particles (hence m SPH = 1O _5 M ), and we set 
Psink — 10 -11 gcm -3 and ^ SINK = 2. The initial density 
profi le is obtained by integrating t he Isothermal Equation 
(e.g. IChandrasekhar fc Wareslll949h from f = to f = 6.45, 
and scaling the result to give total mass M and radius 



10 3 AU. The initial conditions are then generated by relax- 
ing a periodic cube of SPH particles to produce a uniform 
glass, cutting a sphere containing 10 5 particles from this 
cube, and stretching the particle positions to produce this 
density profile. The particles start from rest. 

The cloud collapses to form a single sink surrounded by 
a small accretion disc. Fig. [2] shows, for all the SPH particles 
within ~ 20 AU of the point-mass, the tangential viscous ac- 
celeration, a TV , the radial hydrodynamic acceleration, a RH , 
and the density, p, all as functions of distance from the point- 
mass, shortly after sink creation. The upper row gives results 
obtained with a NewSink, showing that the SPH particles 
near R s (=2AU) have smoothly varying a TV , a RH and p, 
due to the presence of SPH particles in the interaction-zone; 
in particular, a RH is positive (i.e. outward, therefore resisting 
accretion, as it should). The lower row gives results obtained 
with an OldSink, showing that SPH particles near the sink 
boundary experience artificially reduced a TV , a RH and p, due 
to a lack of neighbours in the exclusion-zone; all three fac- 
tors act to increase the net inward acceleration, and thereby 
to artificially increase the rate of accretion. In addition, be- 
cause the SPH particles accreted by an OldSink carry their 
angular momentum with them, rather than transferring it 
to the SPH particles a bit further out, these SPH particles 
can more rapidly move inward, to be accreted in their turn. 
When this test is performed with UrSinks, the results are 
chaotic, and involve the creation of many different sinks; we 
do not discuss these results. 

Fig. [3] shows the evolution of the sink mass, M s , and 
the sink accretion rate, M s , as functions of the time elapsed 
since sink formation, showing that the artificially increased 
accretion obtained with an OldSink is not a transient effect. 
By the end of the simulation, which corresponds to ~ 100 or- 
bital periods at the sink boundary, the OldSink has M s c± 
0.37 M , whereas the NewSink only has M s ~ 0.24 M . 
Thus a numerical simulation using OldSinks overestimates 
accretion rates. The effect on final protostellar masses is 
hard to call. An enhanced accretion rate would mean en- 
hanced radiative and mechanical feedback (if these were in- 
cluded in the simulation) , and this might actually terminate 
accretion earlier, leading to a lower final mass. 

5.3 Boss-Bodenheimer test 

The Boss-Bodenheimer test ([Boss &: Bodenheimerlll979h is 
a standard t est, applied to m any star formation codes. 
For example, iBate et all (|l995h used it to test the original 
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Radius [AU] Radius [AU] Radius [AU] 

Figure 2. (a,d) Tangential viscous accelerations, (b,e) radial hydrodynamic accelerations, and (c,f) densities, for all the SPH particles 
within ~20 AU of the point-mass, shortly after sink creation in the Rotating Bonnor-Ebert Sphere Test. The results on the top row have 
been obtained with a NewSink, and those on the bottom row with an OldSink. 




Time — 1 [yr] Time - t [yr] 

Figure 3. (a) The sink mass, and (b) the sink accretion rate, in the Rotating Bonnor-Ebert Sphere Test, as a function of time elapsed 
since the formation of the sink at to, using a NewSink (red solid lines) and an OldSink (black dashed lines). 



sink algorithm, although of necessity they could only use a cloud having mass M , radius 2 x 10 3 AU, and density field 
rather small number o f SPH particles (A/" SPH ~ 8 x 10 3 ); 

bate fe Burkertl (jl99/t ) repeated the test with ten times as p(r,0,<j>) = 1.6 x 10~ 17 gcm~ 3 {1 + O.5cos(20)} , (46) 
many SPH particles (JV SPH — 8 x 10 4 ), but did not use sinks. 

where (r, 6, <p) are spherical polar coordinates. It is in solid- 
The version we consider here starts with a spherical body rotation with angular speed Q Q = 1.56 x 10 -12 s -1 , 
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Figure 4. A time-sequence of false-colour images of the Boss-Bodenheimer test, performed with NewSinks, -^ sink = 2, p SINK = 
10 — 11 gem -3 and A/" SPH = 102400. Each frame has dimensions (2.4 x 10 3 AU) 2 , and the time (in Myrs) is given in the top-left corner. 
The colour encodes log 1Q (p/g cm -3 ), where p is the density-weighted mean density along the line of sight, and black dots mark sinks. 



and the ratios of thermal and rotational energy to gravita- 
tional energy are a — 0.25 and f3 — 0.20. The gas obeys a 
barotropic equation of state of the form 



T = 10 K <^ 1 + 



io- 



= gem - 



2/5 



(47) 



corresponding to a molecular gas with ratio of specific heats 
7/5 that is approximately isothermal at low densities and 
approximately adiabatic at higher densities. The minimum 
Jeans mass for this equation of state is M MIN = 0.025 M . 

This test is performed with NewSinks and OldSinks, 
using X SINK = 2 and 4 (see Eqn. []]), p SINK = 10~ 13 , 10~ 12 , 
and 10" n gcm- 3 , and A/" SPH = 6400, 12800, 25600, 51200, 

this 



102400, 204800 and 409600. Putting m SPI? = M /./V SPH 



gives 2.4 x 10 



' M <m 

1V± r^" L S 



< 1.6 x 1O~ 4 M . Using 



R s _ X SINK ?7 



lAUX q 



1/3 



PsiNK 



10 -11 g cm" 



-1/3 



(48) 



N \~ 1/s 



we see that R s ranges from ~ 1.3 AU (small X SINK , large 



A4, 



;e p SINK ) to ~ 46 AU (the opposite extremes) 



Since, in all cases p SINK > 10~ gem - , sinks are only 
created in this test once the gas is well into the adiabatic 
regime. 

The initial conditions are set up by cutting a sphere of 
10 5 SPH particles from a relaxed periodic cube, and scaling 
the mass and radius to M and 2 x 10 3 AU. Next, using 
a spherical polar coordinate system, (r, 0,0), with the pole 
along the z-axis, the azimuthal angle of each SPH particle, 
i, is changed from fa to (/>•, where fa + 0.5sin((/>-) cos(fa) = 
(j>i. Finally, each SPH particle is given an initial velocity 
Vi = ^ Q e 2 XTi. When the test is performed with UrSinks, 
the results are chaotic, and the results show no clear trends. 



Since the Boss-Bodenheimer test is intended to illustrate the 
improved stability of the new procedures for sink evolution, 
the results with UrSinks are not discussed further. 

The evolution is illustrated on Fig. [4] which shows the 
simulation performed with NewSinks, ^ sink = 2, p SINK — 
10 -11 gem -3 and A/" SPH = 102400. As the cloud collapses, it 
flattens (due to rotation) and at the same time the azimuthal 
density perturbation is amplified to produce a symmetric bi- 
nary system with a dense filament between the two compo- 
nents. In Fig. [5] we plot the mass, at time t — 0.030 Myr, of 
the first sink to form, against the number of SPH particles, 
jV SPH , for all the simulations performed. Results obtained 
with NewSinks are given in red, and those obtained with 
OldSinks in black. Results obtained with different values of 
Psink are connected by different line styles. Results obtained 
with ^ SINK = 2 (4) are shown in the upper (lower) panel. 
The second component of the binary system forms almost 
simultaneously with the first, and its mass at t = 0.030 Myr 
is almost identical to that of the first component, so it is 
omitted from these plots, to keep them simple. 

With OldSinks, the results are strongly dependent on 
numerical resolution, m SPH = M /A/" SPH , and on the user- 
defined parameters of sink creation, p SINK and ^ SINK . De- 
creasing p SINK causes sinks to form earlier. Increasing X SINK 
and/or decreasing p SINK and/or decreasing A/" SPH makes 
sinks larger (increases R Sl see Eqn. I49[h and this makes 
the artifical enhancement of the accretion rate greater (see 
below for the explanation). Both effects (earlier formation, 
greater artificial enhancement of the accretion rate) increase 
the mass of the OldSink by 0.030 Myr. 

With NewSinks, there is almost no dependence on 
Psink or ^sink' an d provided the minimum Jeans mass is 
well resolved (A/" SPH > 10 5 ), almost no dependence on A/" SPH . 

The reasons why increasing R s amplifies the artificial 
enhancement of the accretion rate onto an OldSink are 
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Figure 5. The mass of the first sink to form in the Boss- 
Bodenheimer Test, recorded at time £ = 0.03Myr, plotted against 
the number of SPH particles, A/" SPH , for three different sink for- 
mation densities p SINK =10 -13 , 10 -12 , and 10 -11 gcm -3 . Re- 
sults obtained using NewSinks are shown with red diamonds and 
red connecting lines, those obtained using OldSinks are shown 
with black stars and black connecting lines. The upper panel is 
for ^ SINK = 2 (the default option), and the lower panel is for 
X —4 
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threefold. First, to be accreted by an OldSink an SPH par- 
ticle has to enter its exclusion-zone, which is evidently eas- 
ier if the exclusion-zone is larger. Second, the flow of SPH 
particles into the exclusion-zone is amplified artificially by 
the steep non-physical gradients at R s . Third, once an SPH 
particle is accreted, it takes with it specific angular momen- 
tum of order (GMsRs) 1 ^ 2 . The larger R s is, the larger the 
amount of angular momentum that is removed by accretion, 
when this angular momentum should be transferred to the 
matter that is next in line to be accreted, thereby reducing 
its chances of being accreted in the immediate future. 

With NewSinks, this does not happen. First, the flow 
of SPH particles into the interaction-zone is not amplified ar- 
tificially by steep non-physical gradients at R s . Second, the 
SPH particles that are accreted by the point-mass have nor- 
mally transferred much more angular momentum (to other 
SPH particles) in the process of migrating into the centre of 
the interaction-zone. Third, whatever angular momentum 
they do add to the point-mass is then returned to the SPH 
particles remaining in the interaction-zone, so that the an- 
gular momentum of the point-mass remains small. 



5.4 Turbulent core collapse 

In order to establish how well NewSinks perform un- 
der more realistic conditions, we simulate the collapse of 
a turbulent prest ellar core, using initial conditions from 
IWalch et all (|20ll specifically core 'A-MM8' with seed 200). 
The core has mass 1.28 M , initial radius 2.5 x 10 3 AU, and 
the density prodile of a critical Bonnor-Ebert Sphere (i.e. 
one with dimensionless boundary radius £ B = 6.45). The ini- 
tial turbulent velocity field subscr ibes to a power spectrum 
of the form P(k) oc k~ 4 (see also IWalch etaflbuld I20l3 ), 
with 2 ^ k ^ 10 (where k—1 corresponds to the core radius). 
The gas is initially isothermal at 11 K, and the equation 
of sta te is evolved with the algorithm of IStamatellos et al.l 
<|2007h , which captures the transport of cooling radiation 
and the dependence of the opacity on density and temper- 
ature. The Mach Number of the turbulence is 1.5, and the 
ratios of thermal and turbulent energy to gravitational en- 
ergy are a = 0.018 and 7 = 0.040. 

The initial conditions are set up by cutting a sphere 
of 157,000 SPH particles from a relaxed periodic cube, and 
scaling the mass and radius so that the inner 128,000 SPH 
particles have total mass M CORE = 1.28 M and outer ra- 
dius R COKE = 2-5 x 10 3 AU; thus all SPH particles have mass 
m SPH = 10~ 5 M . Next the inner 128,000 SPH particles are 
stretched to fit the density profile of a critical Bonnor-Ebert 
Sphere (as in the rotating Bonnor-Ebert Sphere test of Sec- 
tion [5T2J , and allocated a temperature T= 11 K. At the same 
time, the outer 29,000 SPH particles are stretched so that 
they have uniform number-density ten times smaller than 
the number-density of the SPH particles just inside R conE , 
and allocated a temperature T = 110 K. Finally, the tur- 
bulent velocity field is evaluated on a Cartesian grid, and 
the velocities of individual SPH particles are obtained by 
interpolation on this grid. 

The evolution of the core is simulated using NewSinks, 
OldSinks and UrSinks, with three different values of 
Psink = i0 _1 \ 10~ 10 and 10 -9 g cm -3 ; with these high val- 
ues of p SINK the gas is well into the adiabatic heating regime 
by the time sinks are created. In all cases we set AT SINK =2, 
and all nine simulations are terminated at t FIN = 150 kyr. 
Fig. [6] shows the density field on the z — plane at t FIN , and 
the projected positions of the sinks. Table [3] gives the times 
at which the sinks form and their masses at t FIN , plus the 
total mass and number of sinks at t Fm . 

The simulations performed with NewSinks appear to 
be well converged, (i) For all three values of p SINK , three 
sinks are formed, (ii) All three sinks are formed at the same 
time, apart from a small systematic shift that is attributable 
to the fact that as p SINK increases, it takes a little longer for 
the protostar to contract to p SINK . (iii) All three sinks are 
formed in the same place, (iv) The masses of all three sinks 
at t FIN are the same, apart from a small systematic shift 
which is attributable to the fact that as p SINK increases the 
size of the sink decreases, and therefore more of the mass 
destined for the protostar is still in the form of active SPH 
particles. The positions of the sinks and the distribution of 
the residual gas are also essentially the same, again, modulo 
the fact that there is more residual gas that has not yet been 
accreted, when p SINK = 10 _9 gcm -3 , than when p SINK = 
10 _11 gcm" 3 . 

The simulations performed with OldSinks are not con- 
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Figure 6. False-colour images of log 10 (p/g cm -3 ) at time t FIN =150 kyr for the simulations of a collapsing turbulent prestellar core; p 
is the density- weighted mean density along the line of sight, and black dots mark sinks. The upper row shows frames from simulations 
performed with NewSinks, middle row OldSinks, and bottom row UrSinks. From left to right the three columns show results obtained 
with p SINK = 1CT 11 , 1CT 10 , and lO^gcm" 3 . 



verged. Only the times of formation of the first two sinks are 
more-or-less independent of p SINK . For p SINK = 10 -9 g cm -3 
the third sink does not form at all (at least, not by 150 kyr), 
and in all cases the masses at 150 kyr do not vary with p SINK 
in a systematic way. 

With UrSinks , the situation is chaotic, in the sense 
that - apart from the time of formation of the first sink, 
which is independent of the sink type - there is no pattern 
to the changes in creation time and final mass produced by 
changing p SINK . The only systematic trend is that the total 
number of sinks formed increases with decreasing p SINK . The 
positions of the sinks, and the distribution of the residual 



gas are also very different for different values of p SINK - In 
short, the simulations with UrSinks are divergent. 



6 DISCUSSION 

All three types of sink use the density threshold and over- 
lap criteria (Eqns. [2]& [3} for sink creation. NewSinks and 
OldSinks additionally use the potential minimum and Hill 
Sphere criteria (Eqns.|4]&[5]), whereas UrSinks use the ac- 
celeration divergence, velocity divergence, and non-thermal 
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(gem- 3 ) 10- 11 10- 10 10- 9 10- 11 10- 10 10- 9 10- 11 10- 10 10- 9 



ID i — NewSinks — > < — OldSinks — > < — UrSinks — > 

1 86.8,0.12 87.5,0.11 89.5,0.10 86.8,0.16 87.5,0.14 89.5,0.19 86.8,0.13 87.5,0.20 89.5,0.20 

2 87.9,0.12 89.2,0.11 91.2,0.10 88.1,0.14 89.3,0.22 91.2,0.18 87.9,0.11 89.2,0.13 113,0.18 

3 110,0.19 112,0.16 113,0.14 110,0.18 112,0.12 110,0.08 112,0.11 

4 113,0.07 125,0.04 

5 114,0.05 

6 114,0.03 

7 116,0.07 

8 116,0.05 

9 131,0.03 



M TOT 0.43 0.38 0.34 0.48 0.48 0.37 0.62 0.48 0.38 

AL nrr 333 332 942 



Table 3. Properties of the sinks formed in the nine simulations of a collapsing turbulent prestellar core that are depicted in Fig. [6] Each 
column represents a different simulation (different combination of sink algorithm and/or p SINK ), and lists sequentially the formation 
time, in kyrs, and the mass, in M , at t FIN = 150 kyr (thus, for example, in the simulation using NewSinks and p SINK =10 -10 gcm -3 
the third sink is created after 112 kyr and by 150 kyr has M s =0.16 M ). In addition, we give, at the foot of the table, the total mass of 
sinks formed, M TOT , by £ FIN =150 kyr, again in M , and the total number of sinks, A/" TOT . 



energy criteria (Eqns. [33] [34] & [35]). As a consequence, 
UrSinks are prone to artificial sink creation. 

A NewSink comprises a point-mass and a concentric 
spherical interaction-zone populated by SPH particles; these 
SPH particles are accreted by the point mass on an extended 
timescale (regulated accretion) and the angular momentum 
that they add to the point-mass is returned to the SPH parti- 
cles remaining in the interaction zone, again on an extended 
timescale. In contrast, an OldSink or UrSink comprises a 
point-mass and a concentric spherical exclusion-zone that 
usually is almost devoid of SPH particles. Most SPH parti- 
cles entering the exclusion-zone are immediately assimilated 
by the point-mass (instantaneous accretion), taking their an- 
gular momentum with them, and this results in artificially 
steep gradients in the vicinity of the exclusion-zone bound- 
ary, and artificial removal of angular momentum. 

The potent i al mi nimum criterion was introduced by 
iFederrath et al.1 ([20101 ), and is effective because the gravi- 
tational potential derives from all the matter in the compu- 
tational domain, and therefore is not strongly influenced by 
particle noise. In combination with the Hill Sphere criterion, 
it stops density peaks that are destined to be sheared apart 
by some external gravitational field from forming sinks. 

Regulated accretion is designed to cure the problems 
that arise from the lack of SPH particles in the exclusion 
zone of a standard sink (OldSink or UrSink) and the non- 
physical flow of angular momentum into a standard sink. 

The tests we have implemented here are not intended to 
be exhaustive, but to give an indication of the improvement 
in performance that NewSinks deliver. First, NewSinks 
ensure relatively well-behaved hydrodynamics at the sink 
boundary. Second, they do not act as sinks of angular mo- 
mentum. Third, as a consequence of these first two features, 
they predict approximately converged results (i.e. creation 
times, masses, accretion rates, etc.) that are more-or-less 
independent of the user-defined sink parameters. 

We stress that convergence can only be tested using 
structures that develop from prescribed perturbations, i.e. 
perturbations that can be accurately reproduced in the ini- 



tial conditions of simulations having different resolutions 
and/or different sink parameters. Thus, for example, in the 
Boss-Bodenheimer Test, when the two original protostars 
and their accretion discs interact violently at periastron 
(~ 0.032 Myr) several new sinks may form. However, these 
have their origin in fluctuations due to particle noise, and 
therefore their genesis is a chaotic process that cannot be 
used as the basis of a convergence test. It is for this reason 
that we chose to measure the mass of the first protostar to 
form, and to make this measurement at ~ 0.030 Myr. 

NewSinks incur a significant computational overhead, 
but, given the non- convergence and non-physicality inherent 
in the use of standard sinks, this overhead is a price that 
has to be paid. For example, in the turbulent core collapse 
test of Section switching from OldSinks to NewSinks 
increases the number of CPU-hours from 68 to 253 when 
10 -11 gcm -3 , and from 313 to 766 when p SINK = 
3 . The increase in CPU time derives mainly from 
the fact that the SPH particles inside the interaction-zone 
are - of necessity - evolved with very short timesteps. 

A further advantage of NewSinks is that, since they 
produce converged - and therefore hopefully accurate - ac- 
cretion rates, one can have more confidence in using these to 
generate accretion luminosities (radiative feedback) and out- 
flows (mechanical feedback) from protostars or black holes. 
Sub-grid physics can easily be added to a NewSink, for 
example a model of stellar evolution with accretion (cf. 
Krumholz et al.l 120071), a model for episodic accretion (cf. 



PsiNK 

10~ 9 gem 



Stamatel 



os et al.ll2011 ), or a model with mechanical feed - 



back (cf. IStamatellos et al .11 20051 : ICunningham et alJlioITh . 

The NewSi nk algorithm has som e similarities to the 
scheme used by ISpringel et al.l (|2005h to model accretion 
onto black holes in galactic nuclei. This scheme estimates 
an accretion rate using the Bondi-Hoyle-Littleton theory 
(iHovle fe Lvttletonlll939l : iBondi fe Hovlelll944l : lBondilll952h 
and the local density, sound speed and velocity dispersion; 
the SPH particles to be accreted are chosen stochastically, 
but with a weighting that preferentially selects those near- 
est the BH; and the dynamical mass of the BH is increased 
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smoothly so that the accretion of small numbers of massive 
SPH particles does not cause large feedback fluctuations. We 
have also experimented with using an average of the Bondi 
accretion rates for the SPH particles in the interaction-zone 
to determine the time-scale for radial accretion, but we find 
that Eqn. @ is better behaved. Moreover, it may some- 
times be more realistic to treat accretion onto a BH in a 
galactic nucleus using the disc-accretion limit (as would hap- 
pen automatically w ith a NewSink). Additionally, in the 
ISpringel et al.l (|2005h scheme the BH assimilates all the an- 
gular momentum of an accreted SPH particle, which may 
lead to unrealistic effects. 

We stress that the issues of robust sink creation that we 
have addressed here are quite distinct from those of resolu- 
tion, which we have not addressed. For example, if, as in our 
code, inter-particle gravity is kernel-softened, the minimum 
Jeans mass is only adequately resolved if it contains very 
many more particles than the mean number of neighbours 
(i.e. a large multiple of A/" NEIB ), so that gravitational inter- 
actions between the majority of particles are not softenned; 
that is a resolution issue. However, if sinks are not created 
and evolved in such a way as to avoid spurious creation, 
unphysical boundary gradients, and unphysical assimilation 
of angular momentum, then, even if the self- gravitating gas 
dynamics that leads to sink creation and subsequently feeds 
sinks is faithfully simulated, with high resolution, the prop- 
erties of the sinks formed are not converged, and cannot be 
trusted; these are the issues that NewSinks deal with. 



7 CONCLUSIONS 

We have developed a new algorithm for the creation and evo- 
lution of sink particles in SPH, which has four user-defined 
free parameters, (p SINK , 

^sinio *hill» <*ss)» with default 
values (10 -11 gcm -3 , 2, 4, 0.01). This algorithm has sev- 
eral advantages over the standard algorithms used hereto- 
fore in SPH simulations of star formation, (i) SPH particles 
flow into the volume of a NewSink (its interaction-zone), 
and are only accreted by the point-mass later, piecemeal, 
and usually over many timesteps. Consequently the hydro- 
dynamics in the vicinity of the boundary of a NewSink 
is not so severely corrupted by poorly evaluated gradients, 
(ii) The angular momentum of SPH particles assimilated by 
the point-mass is transferred back to the SPH particles in 
the interaction-zone, so that the point-mass does not act 
as a sink for angular momentum; this angular momentum 
is in any case much reduced by the fact that the accreted 
SPH particles have normally had to lose angular momen- 
tum to other SPH particles before they become the closest 
SPH particle to the point-mass (and therefore next in line 
to be accreted). These factors mean that the rate of accre- 
tion onto a NewSink is both better informed by physics 
than that onto a standard sink, and relatively independent 
of the user-defined parameters of sink creation. Provided the 
simulation has sufficient resolution, well converged results 
are obtained up to density p SINK - We have measured the 
new algorithm against four tests: Bondi accretion, the col- 
lapse of a uniformly rotating Bonnor-Ebert sphere, the Boss- 
Bodenheimer Test, and the collapse of a turbulent prestel- 
lar core. In all cases the NewSinks appear to produce quite 
reliable results (i.e. accurate where there is a known solu- 



tion, and otherwise approximately converged), whereas the 
standard algorithms do not. If p SINK is fixed, then a sim- 
ulation with NewSinks requires a factor between two and 
four times more computing resource than one with stan- 
dard sinks,. However, this expense must be set against the 
fact that results obtained with NewSinks are more credible. 
Furthermore, since the results obtained with NewSinks are 
converged at lower values of p SINK , one can safely reduce the 
computing time - by a comparable factor - by decreasing 
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